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We discuss the reliability of integral-equation methods based on several commonly used closure 
relations in determining the phase diagram of coarse-grained models of soft-matter systems char¬ 
acterized by mutually interacting soft and hard-core particles. Specifically, we consider a set of 
potentials appropriate to describe a system of hard-sphere colloids and linear homopolymers in 
good solvent, and investigate the behavior when the soft particles are smaller than the colloids, 
which is the regime of validity of the coarse-grained models. Using computer-simulation results 
as a benchmark, we find that the hypernetted-chain approximation provides accurate estimates of 
thermodynamics and structure in the colloid-gas phase in which the density of colloids is small. On 
the other hand, all closures considered appear to be unable to describe the behavior of the mix¬ 
ture in the colloid-liquid phase, as they cease to converge at polymer densities significantly smaller 
than those at the binodal. As a consequence, integral equations appear to be unable to predict a 
quantitatively correct phase diagram. 


I. INTRODUCTION 

Integral-equation methods are a very powerful tool to determine the thermodynamics and the liquid 
structure of simple fluids 0 , 0 . They rely on different approximate closure relations which, supple¬ 
mented by the Ornstein-Zernike (OZ) equation, allow a direct and numerically fast determination of the 
pair correlation functions as well as of thermodynamic quantities like pressure, compressibility, chemical 
potential, ... For simple fluids these methods cannot compete nowadays with Monte Carlo and molecular- 
dynamics simulations. Nonetheless, they have the advantage of providing reasonably accurate estimates 
of thermodynamic quantities with a very limited effort, and they are therefore a very valuable tool when 
the system under investigation depends on many parameters, for instance in the case of multicomponent 
systems. Moreover, they are still very useful for the analysis of systems for which atomistic simulations 
are particularly slow, for instance in glassy systems; see, e.g.. Refs. 01. 

Liquid-state integral equations have also been extensively used to compute fluid-fluid phase-coexistence 
lines. In the density region in which the system demixes, integral equations may not converge, or may 
converge to physically unacceptable solutions. The relation between the boundary of this nonconvergence 
region (we will call it termination line) and the binodal and the spinodal curves characterizing the two- 
phase unstable region has been the subject of many studies, see, e.g.. Refs. il. In particular, it has 
been shown that, except in the case of very simple approximations, thermodynamical quantities do not 
show any particular divergence on this line, hence it cannot be taken as an approximate estimate of the 
spinodal line. However, it is usually assumed that it is somewhat close to the line where phase separation 
occurs. 

In this paper we wish to investigate the reliability of integral-equation methods for the determination 
of the phase diagrams of typical coarse-grained models of soft-matter systems. We consider here a binary 
mixture of soft and hard spheres of different sizes with an intrinsic nonadditive nature. Although we 
take specific pair potentials, appropriate to describe, in a coarse-grained fashion, a binary system of 
hard-sphere colloids and long polymers under good-solvent conditions [l3, [Hi, the conclusions should 
apply to a general class of soft-matter systems that can be modelled as mixtures of soft and/or hard 
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spheres, interacting via short-range potentials (l2l - [l^ . The phase diagram of the coarse-grained model 
has been accurately determined in Ref. [l^, by means of Monte Carlo simulations, for different values 
of the polymer-to-colloid size ratio. Here, we investigate the same problem by using integral-equation 
methods. We employ the hypernetted-chain (HNC) , the Percus-Yevick (PY), the Rogers-Young (RY), 
and the reference HNC (RHNC) closures [ll. fisl - l^ . For each of them we determine the termination line, 
whose position is then compared with the Monte Carlo binodal with the purpose of understanding if this 
line provides a reasonable approximation of the boundary of the two-phase region. For small polymer 
densities, we will also be able to compute by Monte Carlo simulations the bridge functions—quantities 
that have an intrinsic interest in liquid-state theories—which can then be compared with the approximate 
ones considered in the different approaches. 

The paper is organized as follows. In Sec.|ll]we define the model, report the definitions of the different 
closures we use, and the explicit expressions of the quantities that are considered in the paper. In Sec. im 
we present our results. In Sec. IHI Al we determine the termination line for the different closures for two 
different values of the polymer-to-colloid size ratio q, q = 0.5 and q = 0.8. In Sec. IHI HI we compare the 
integral-equation predictions for structure and thermodynamics with Monte Carlo results. In Sec. IHI Cl 
we determine the bridge functions with Monte Carlo methods and compare them with those used in the 
different integral-equation approaches. In Sec. IHI PI we consider a novel approximation that uses the 
Monte-Carlo determined bridge functions. Finally, in Sec. IIVI we draw our conclusions. Technical details 
are reported in Appendix The explicit expressions of the potentials are reported in Appendix |B] 


II. DEFINITIONS 


A. The model 

We consider a mixture of mutually interacting hard spheres of radius Rc and of soft particles with a 
typical interaction range Rg. Specifically, we consider here a set of potentials which are appropriate to 
describe a system of hard-sphere colloids of size Rc and linear homopolymers in good solvent of radius of 
gyration Rg, after tracing out the monomer degrees of freedom and replacing each chain with a particle 
coinciding with its center of mass. The coarse-grained model is accurate only if polymers are dilute, i.e., 
for = iTrNpRg/{3V) < 1 {Np is the number of colloids in the volume V), and if the polymer-to-colloid 
size ratio q = Rg/Rc satisfies q <1 discussion of the accuracy of the model, with a comparison 

with full-monomer results is presented in Ref. IS)- 

Polymer-colloid solutions have been extensively studied [2ll - [^ . because of their rich phase diagram, 
which presents fluid-fluid and fluid-solid coexistence lines, and because of their technological relevance 
[^. In this paper, we will not be interested in using the model to predict their phase behavior. Rather, 
we take it as a typical soft-matter system and use it as reference model for which we can study the 
predictivity of the different closures that are typically used in integral-equation studies. We will consider 
three different values of g = ^/Rc, q = 0.5,0.8, and 1. The corresponding pair potentials have been 
determined in several papers Here we shall use the accurate scaling-limit results of Refs. [^1^. 




FIG. 1: Left: polymer-polymer pair potential pVpp{b) as a function of 6 = r/Rg\ right: polymer-colloid potential 
PVcp{b;q) as a function of 6 = r/Rg for q = 0.5, 0.8, 1. For q — 0.8 and q = 0.5, we assume Vcp{b;q) = oo for 
b < 0.90, 1.91, respectively. 
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They are reported for completeness in Appendix [B] and shown in Fig.[TJ The polymer-polymer potential 
is essentially Gaussian with Vpp{b = 0) ~ l.SfcsT at full overlap. On the other hand, the nature of the 
polymer-colloid potential depends on q. For q < 1 the potential is expected to be infinite at full overlap 
and very large for r < Rc = Rg/q- Then, it decays fast, with a tail that is small for r > Rc + 2Rg. Note 
that we have not been able to determine liVcp{b\q) in the small-r region in which j3Vcp{b]q) > 10. For 
these values of r we simply assume l3Vcp{b; q) = -|-oo for q = 0.5 and q = 0.8. For q = 1 we performed a 
linear extrapolation. 


B. Closure relations 


In the integral-equation approach the basic ingredients are the pair correlation functions hap (r) (a and 
P label the two species) and the direct correlation functions Cap(r). They are related by the Ornstein- 
Zernike (OZ) [ij relations 

hajuih^ — T E Ca7 ik)p^h^p{k), (1) 
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where we denote with f{k) the (three-dimensional) Fourier transform of any function /(r). To compute 
the quantities of interest, the OZ relation must be supplemented by a closure relation, which can be 
written in general form as 

9a0{r) = exp[hap{r) - Cap(r) + bap{r)\, (2) 


where Vap{r) are the pair potentials, ga/sir) = hap{r) -I- 1 is the pair distribution function and bap(r) is 
the so-called bridge function. The latter quantity cannot be computed exactly, hence we consider several 
different approximations: 

(i) Hypernetted chain (HNC) closure [I|. We simply set bap{r) = 0 for all a,p. This approximation is 
very accurate for soft potentials [T|. 

(ii) Mixed HNC/Percus-Yevick (PY) closure. For hard spheres the PY closure relation [I| 

gMr) = + hapir) - Capir)] (3) 


is more accurate than the HNC closure. Here we consider the HNC closure for polymer-polymer 
and polymer-colloid correlations and the PY closure for the colloid-colloid correlations. 

(hi) Rogers-Young (RY) closure [ll,|33. This closure mixes the HNC and the PY closures, adding free 
parameters that are tuned to obtain thermodynamic consistency. It is defined by 


gapir) 




, exp[{hap{r) - Cap{r))fap{r)] - 1 

U{r) 


( 4 ) 


where the function fapir) is given by 


/a/3 




( 5 ) 


Note that, for Xap —!> 0 we recover the PY closure, while in the opposite limit, Xa/3 oo, we 
reobtain the HNC closure. In most of the discussion we have considered a single optimization 
parameter, setting Xap = x/sa/3, Scc = Ro Spp = Rg, Spc = {Rc + Rg)/‘2- The parameter x has 
been determined as discussed below. 


(iv) Reference HNC (RHNC) closure In this approach one sets bap{r) = b^p {r\ Rp, Rc), where 

the latter quantities are the bridge functions of a system of additive hard spheres of radii Rp and 
Rc at the same densities of the polymers and colloids in the original system. The polymer effective 
radius Rp is determined by using the Lado criterion 


y^^XgXp 

0.(5 



r^dr [hap{r) - hg^ (r; Rp, Rc)] 


db(^^{r-, Rp, Rc) 
dRp 


= 0 , 


( 6 ) 


where Xg = Na/{Nc + 
as discussed in Refs. |T 


Pa/{Pp + Pc)- The bridge functions b((g {r; Rp, Rc) can be computed 


Solving simultaneously the OZ and the closure relations, one obtains hap{r) and Cap{r). Then, one can 
use them to compute thermodynamic quantities. 
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C. Observables 


We will be interested in computing the pressure. One possibility consists in using the virial expression: 


= p 



(7) 


where /3 = l/fc^T, p = pp + pc, and the quantities Zap are given by 


^aP — 


2tt 

~l^PaPP 

3p 



dPVap 

dr 


9ap{r). 


( 8 ) 


This expression cannot be applied to hard spheres, since the potential is discontinuous. In this case we 
have 

= -^plgac{2Rc). 9 

3p 

Eq. ([5]) is also not convenient in the polymer-colloid case as pVpc{r) diverges as r —0. In the HNC case, 
this problem can be overcome by rewriting Eq. m as 


Zap — 


r'^dr ■ 


de 

dr 


,hci3(r)-c^l3{r) 


( 10 ) 


A similar formula can be analogously obtained in the case of the RY closure. 

Another quantity we shall be interested in is the isothermal compressibility kt that can be either 
computed by using the virial route 



“ V dpp dpc )p/"" 



'y ^ PaPpCap{0)- 
0.(5 


( 11 ) 


( 12 ) 


The two expressions are thermodynamically equivalent. However, when an approximate closure is used, 
two different results are obtained, as a consequence of the thermodynamic inconsistency of the approach. 
In the RY case, the parameter y is fixed so that the two different routes provide the same result for kt- 
Einally, we shall consider the structure factors 

Sapi.k') — ^a/3 4“ PaPphapi.^) ^ (18) 

and the concentration structure factor 


Scik) = XpXc [xpSccik) + XcSppik) - 2^XpXcScpik)] . (14) 

For fc —>• 0, 1/Sc{k) —>■ d'^l3g{xp,P)/dXp, where g{Xp,P) is the Gibbs free energy per particle. Hence, its 
divergence signals the thermodynamic instability of the homogeneous phase. 


III. RESULTS 

In order to solve the coupled integral equations, the correlation functions are discretized on a regular 
grid. We usually take a step size Ar/Rg = 10“^ and truncate the correlation functions at Rmax/Rg = 
NAr, with N = 32768. As we discuss in appendix lAl these choices make truncation and discretization 
errors negligible. We use the standard Picard iterative method, which converges quite fast, except close 
to the termination line. We improve convergence by considering a mixing parameter a. If c|^^(r) and 
Cg"j(r) indicate the direct correlation functions at the beginning and at the end of the n-th step of the 
iterative procedure, respectively, we set c|^~'’^^(r) = (1 — Q;)cj(((^(r) -1- ac^^^ir). Far from the termination 
line, a is not a relevant parameter. However, close to the termination line, convergence is only obtained 
if a is small. In some cases, we took a ^ 10“^. 
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FIG. 2: Results for q = 0.8 and the HNC/PY closure (in this case there is no termination line). Left: isobars 
corresponding to = 1, 2, 3, 4, and 5 in the 'he, 4>p plane. We also report the Monte Carlo estimates of 

the binodal (empty squares, MC). Note that isobars with /SPR^ > 3 go through the two-phase region. Right: 
1/Sc{k = 0) along the same isobars reported in the left panel as a function of Xp = Np/{Np + Nc). 


A. Termination lines 

In order to determine the termination line, we work as follows. We fix the colloid volume fraction 
($c = AttRINc/^V, where W is the number of colloids present in the box of volume V) and solve the 
equations for a small value of the polymer density. Typically, if <f>p = 4:TrRgNp/3V {Np is the number 
of polymers present in the box of volume V), we start at $p si 0.005 for ^ 0.2 and at <l>p ~ 0.01 for 
smaller colloid volume fractions. Then, we increase by steps A$p. 

For g = 1 we have been able to increase up to 2.5 for all values of $c < 0.45: We always find a 
regular solution of the integral equations. This is not surprising, as, for this value of g, Monte Carlo 
simulations indicate that the fluid-fluid binodal either does not exist or is located at quite large values 
of the polymer volume fraction. In particular. Ref. [13 found no phase transition up to = 2.12, 1.73, 
1.33 for = 0.1, 0.2, 0.3, respectively. 

For g = 0.8 integral equations do not show any singular behavior if the HNC/PY closure is used. Also 
in this case we have been able to solve the equations for any <l>p < 2.5 and 4>c < 0.45. No phase demixing 
is observed, as it is evident from the behavior of 1/Sc{k — 0) along five different isobars shown in Fig. [5] 
At the critical point 1/Sc{k = 0) should vanish. Instead, it increases as the pressure P is increased, 
with no indication of a zero for some values of P and Xp. Apparently, the HNC/PY closure fails even in 
reproducing the qualitative behavior of the system. 


TABLE I: For each q and 'Fc (first two columns), we report the polymer volume fraction at which integral 
equations no longer converge to a physical solution for three different closures: HNC, HNC/PY, and RY. In the 
last column we report the polymer volume fraction at which the binodal, as computed by Monte Carlo 
simulations [T^], occurs. We have also computed the termination line for the RHNC closure, for q — 0.5 and 
•Fc = 0.3: -Fp = 0.104. 


q c&c 

HNC HNC/PY 

RY 

^bin 

0.5 0.10 

0.87 

> 2.5 

0.88 

0.69 

0.20 

0.23 

> 2.5 

0.34 

0.53 

0.30 

0.090 

0.15 

0.18 

0.38 

0.40 

0.036 

0.07 

0.115 

0.255 

0.8 0.10 

> 2.5 

> 2.5 

> 2 


0.20 

0.61 

> 2.5 

> 2 

1.0 

0.30 

0.175 

> 2.5 

0.39 

0.75 

0.40 

0.086 

> 2.5 

0.29 

0.5 < < 0.6 


For g = 0.8 a termination line is observed if we use the HNC or the RY closures, while for g = 0.5 
a no-convergence domain is observed also by using the HNC/PY closure. Results for the termination 
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FIG. 3: Estimate of Sc{k) for >l>c = 0.3 and <l?p = 0.0905, on the termination line. Here q — 0.5 and we use the 
HNC closure. 
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FIG. 4: Estimate of Sc{k — 0) for three different A4>p and the HNC closure. We start from the values of Cap{r) 
and ha:p{r) for = 0.9 and increase <E>p by steps A<l?p up to iHp = 0.0915, then we decrease <E>p with the same 
schedule up to = 0.086. The termination line occurs for = 0.0903. 

lines for both values of q are reported in Table ID The termination line is determined as follows. Starting 
from the initial value we subsequently solve the equations for <l>p”^ = + nA^p, starting the 

iterations for the n-th density from the solution at $p” If A<i>p is large or the mixing parameter in 
the Picard iterations is of order 1, we end up at a density $p^^ where the iterations no longer converge. 
Then, we consider again the solution at $p^ but now we significantly decrease A$p and the mixing 
parameter (typically we take a parameter as small as 0.01). If we increase again $p, we now observe that 
the Picard iterations always converge. However, at a very specific value of $p the stable solution is no 
longer physical, as Sc{k) becomes discontinuous at a finite value of k. We identify the termination line as 
the smallest polymer density at which S'c(fc) (the same occurs for all structure factors Sapik)) develops a 
discontinuity. An example is shown in Fig.[3l where we report Sc{k) for q = 0.5, = 0.3, $p = 0.0903, as 

obtained by using the HNC closure. It is interesting to observe that while the position of the termination 
line is independent of the protocol used to increment $p, the singular solution depends on A$p. For 
instance, in Fig. |4] we show the estimates of Sc{k = 0) as a function of d)p for three different values of 
A<I>p. Incrementing il>p, at the termination line $p = 0.0903 we always observe a jump in 5'c(0) to a 
new value. However, such value depends on A^bp. If we further increase <I>p beyond the termination line 
and then decrease again il>p, the unphysical solution appears to be stable: The structure factor changes 
smoothly with <I>p. Moreover, once $p is again below the termination-line value, if we use a small mixing 
parameter, we always obtain the unphysical solution. 

The termination lines for the different closures are reported in Fig. [51 In general, we find that the 
RY closure performs better than the HNC one, which stops converging at very small values of $p in the 
colloid-liquid phase. In all cases, however, the termination line is significantly below the correct binodal, 
especially in the colloid-liquid phase ^ 0.25. Clearly, the convergence to an unphysical solution is not 
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FIG. 5: Phase diagram for q — 0.5 (left) and q = 0.8 (right). We report the binodals obtained by Monte Carlo 
simulations (MC), the corresponding critical point (MC-CP), and the termination lines for each of the closures. 


directly related to singularities in the thermodynamic behavior of the model. Therefore, the termination 
line provides a very poor approximation of the phase-separation line. 

Let us finally consider the RHNC closure. Since this approach is quite complex, we have only analyzed 
one case: q — 0.5 and $c = 0.3. For $p « 0, the effective radius Rp is equal to 0.837i?g. This is a 
completely reasonable value, indicating that polymers are effectively equivalent to hard spheres of radius 
approximately equal to Rg. As <i>p increases, the effective radius Rp decreases quite rapidly: for $p = 0.1 
we find Rp = 0.60Rg. Again, this is consistent with intuition, as we expect the polymer to shrink as $p 
increases. Unfortunately, we are not able to go much beyond $p = 0.1, as the RHNC equations cease to 
converge at $p = 0.104. Hence, this approach represents only a modest improvement with respect to the 
HNC approach (the HNC termination line occurs at <l>p = 0.090). 


B. Structural behavior in the homogeneous phase 

We wish now to compare the integral-equation predictions with the Monte Carlo ones in the homoge¬ 
neous phase. We consider the case q = 0.5, in which a termination line occurs for all considered closure 
relations. We begin by analyzing the structure factors Sap{k = 0), which are directly related to thermo¬ 
dynamics by the compressibility equations [H, 13. In Fig. [6] we report the corresponding estimates for 
two values of $c, = 0.1 and 0.3, that lie on opposite sides with respect to the critical point located at 

4*c,crit — 0.25, 4^pjCrit — 0 .46, as estimated by Monte Carlo simulations [II- 

For = 0.1 the HNC/PY closure significantly underestimates the structure factors. Clearly, |S'q-/3(0)| 
increases too slowly as $p increases, explaining why convergence is observed at least up to $p = 2.5, see 
Table m The HNC and RY estimates increase faster. The latter are more accurate than the HNC ones 
for small densities, but they significantly underestimate |5'q,,s(0)| close to the binodal, which is located 
at $p ft! 0.70 [l3|- The fact that the RY results are less accurate than the HNC ones near the binodal 
may be surprising, as the RY closure is a generalization of the HNC closure. It simply indicates that the 
requirement of thermodynamic consistency does not necessarily lead to more accurate results. Note that 
both HNC and RY integral equations also converge for some values of <I>p in the metastable region beyond 
the binodal, see Fig. [5) In this domain the structure factors Sap{Q) are quite large [on the binodal, Monte 
Carlo simulations give 5'pp(0) = 11.5(6), 5'cp(0) = —6.7(4), 5'cc(0) = 4.2(2)]. Therefore, even though we 
do not observe an exact divergence of Sctp{k = 0), for this value of we can take the termination line 
as a good estimate of the spinodal. 

For = 0.3 the behavior is quite different and the termination lines occur at values of 4>p significantly 
smaller than that of the binodal. Moreover, integral equations stop converging when the structure factors 
|5'a^(0)| are relatively small, at least if compared with the values they assume on the binodal at = 0.1. 
For instance, the HNC and HNC/PY equations both cease to converge when 5'pp(0) ft 3, while S'pp(O) ft 5 
on the RY termination line. Comparing the integral-equation estimates with the Monte Carlo results, we 
see that the RY closure is here the most accurate, in agreement with previous studies [H, [3 : although 
it fails to converge well before the binodal. As for the RHNC, the estimates of S'cc(O) and S'cp(O) are 
consistent with the RY ones and the Monte Carlo data up to $p ft 0.08. On the other hand, the RHNC 
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FIG. 6: Structure factors Sa^{k — 0) for q — 0.5 and <I>c = 0.1 (left) and >l>c = 0.3 (right). Lines are the results 
obtained by using the HNC, HNC/PY, and RY closures. Symbols are Monte Carlo data. For = 0.3 we also 
include results for the RHNC closure and results obtained by using the zero-polymer-density Monte Carlo bridge 
functions (MC-B), as discussed in Sec. IIIID] 


estimates of <5'pp(0) increase too fast for > 0.04, looking similar to the HNC estimates. Also in this 
case the termination line occurs for S'pp(O) ~ 3. 

As a second test let us compare the pair distribution functions. For = 0.1 and $p = 0.7, i.e. 
on the binodal, see Fig. [71 all closures reasonably reproduce the polymer-polymer distribution function. 
Deviations are instead observed for the polymer-colloid and especially for the colloid-colloid distribution 
function. The largest deviations are observed for the HNC/PY closure. For instance, the colloid-colloid 
correlation is significantly underestimated at contact. While an extrapolation of the Monte Carlo data 
predicts gcc{‘2Rc) ~ 13-14, we estimate gcc{‘2Rc) ~ 4 by using the HNC/PY closure. The RY closure 
performs better, although it is also unable to predict the correct value of gccij’) at contact and slightly 
overestimates gcpif) at the first peak. As for the structure factors, the HNC closure is the most accurate 
one for this value of <i>c, as the HNC curves fall on top of the Monte Carlo data. 

At = 0.3 the behavior is quite different, see Fig. [71 For $p = 0.085, close to the HNC termination 
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FIG. 7: Pair correlation functions gap{r) as a function of fe = f/Rg for q = 0.5 at two different state points: 
<I>c = 0.1, "Fp = 0.70 (left) and <I>c = 0.3, = 0.085 (right). Lines are the results obtained by using the HNC, 
HNC/PY, RY, RHNC closures. Symbols are Monte Carlo data. 


line, HNC results are not accurate, especially for gpp{r), which is significantly overestimated for r < 2Rg. 
The value of gcc(r) at contact is also significantly overestimated. The HNC/PY closure gives results that 
are only marginally better than the HNC ones, while the RY estimates are in full agreement with the 
Monte Carlo data. The RHNC estimates of gcc{r) and gcp{r) are in agreement with the data, but this 
is not the case for gpp{r), which is overestimated for 1 ;< f/Rg ^ 2, the region in which the correlation 
function shows the first peak. At = 0.15 we only have RY data, as integral equations no longer 
converge for the other closures. The results are reported in Fig. [S] Pair correlations gcc(^) and gcp{r) are 
well reproduced, while relatively small deviations are observed for gpp{r). Apparently, RY estimates are 
relatively accurate even close to the corresponding termination line, located at $p = 0.18. 
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FIG. 8: Pair correlation functions gap{r) as a function of 6 = r/Rg for q = 0.5 at <I>c = 0.3, = 0.15. Lines are 

the results obtained by using the RY closure. Symbols are Monte Carlo data. 


C. Bridge functions at zero polymer density 

The failure of integral-equation methods to reproduce the thermodynamics for ^ 0.2 and to provide 
a reasonably accurate estimate of the boundary of the two-phase region clearly indicates that none of 
the closures we used is appropriate for the problem at hand. To understand better the origin of the 
discrepancies, we now compare the bridge functions used in the integral-equation approach with the 
exact estimates obtained numerically, by using the MC results for the pair correlation functions. For this 
purpose we should compute gapir) accurately on large boxes. It turns out that this is feasible only for 
<l>p —^ 0, the case we will study below. 

The input numerical quantities are gcci'i') (we use the accurate expressions that can be obtained as 
discussed in Refs. jlOilly)’ 9cp{f), and gpp{r). To determine the last two quantities, we perform simula¬ 
tions for different values of on systems of linear size L/Rg = 32, 24 for q = 0.5 and 1, and perform 
an extrapolation to <l>p —>■ 0. Then, we determine the direct correlation functions by inverting the OZ 
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FIG. 9: Bridge functions for g = 1 as a function of 6 = r/Rg-. on the left we report bcp{r), on the right bpp{r). 
Top: "he = 0.1; bottom: = 0.3. We report the Monte Carlo estimates (MC) as well as those obtained by using 

the different closures. RY-2 labels the results obtained by using the two-parameter RY closure discussed in the 
text. 


relations, which, for $p —^ 0, simplify to 


Finally, we define 


dec (^) 

Ccp{k) 

Cpp ( /h ) 


^cc(^) 

1 + Pehceik) 

kep {k) pcCcc{k)hcp{k) ^ 
hpp (k) - Pc ^cp (fc)hcp(fc). 


bapir) = In + Capir) - hap{r). 


(15) 


(16) 


We will focus on the polymer-polymer and colloid-polymer functions, as bedr) depends only on the hard- 
sphere fluid, a case that has already been extensively discussed in the literature. Note that jSVepir) is 
large for r < Rc, so that Papir) is not determined accurately for these distances. Hence, we are not able 
to obtain reliable estimates of bcp{r) for r < Rc- 

For the HNC or the HNC/PY closure, we have hcp{r) = bpp{r) = 0. In all other cases, the bridge 
functions are obtained from Eq. (HU), using the correlation functions obtained by means of the different 
closures. For the values of r for which Vcp{r) is large, it is convenient to express 5 cp(?’)e^^'=*’^’'^ in terms 
of hcp{r) — Ccp(r) using the closure relation. This trick allows us to compute the bridge functions bcp(r) 
inside the core region r d Rc, although here they cannot be compared with the Monte Carlo results. In 
this section we do not consider the HNC/PY, as it has the same bridge functions of the HNC closure. 
We will instead discuss the full PY closure, in which Eq. is used for all correlations. 

The bridge functions for = 0.1 and 0.3 are reported in Figs. |9] and [TO] for q = 1 and 0.5, respectively. 
For = 0.1 the bridge functions are tiny, explaining why the HNC closure works reasonably well. The 
PY and RY closures are essentially equivalent. Small deviations are evident for g = 1 and r < 2Rg — but 
in this range data become increasingly less accurate — while for q = 0.5 no deviations are observed in 
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FIG. 10: Bridge function for q = 0.5 as a function of 6 = 'f'/Rg'- on the left we report bcp{r), on the right bpp{r). 
Top: "he = 0.1; bottom: = 0.3. We report the Monte Carlo estimates (MC) as well as those obtained by using 

the different closures. RY-2 labels the results obtained by using the two-parameter RY closure discussed in the 
text. 


the region in which data appear to be reliable. As increases, the bridge functions become increasingly 
negative for small values of r. For q = 1 and <l>c = 0.3, none of the closures appear to be accurate, 
although the RY closure is marginally better, and large deviations are observed for r < 2Rg. For q — 0.5 
the RY closure reproduces well bcp{r) up to r si 2Rg —the region outside the colloid core. On the other 
hand, deviations are clearly observed for bpp(r) when r < Rg. The PY closure is clearly worse, as it 
underestimates both bridge functions for r < 2Rg-3Rg. 

The RY optimization at = 0 uses only the colloid-colloid correlations. Indeed, in this limit the 
consistency condition is 


9/3p(vir) 

dpc 


Pp — 0 


1 PcCcc(O)- 


(17) 


Therefore, one might think that the relatively poor agreement for the polymer-polymer correlations for 
small values of r is related to the fact that the procedure does not take into account polymer properties. 
We have thus considered a two-parameter optimization. We set Xpp = Xi/Rg Xcc = X'ilRc as 
free parameters, while Xpe is, somewhat arbitrarily, set equal to (xi + X 2 )/(.Rg + Rc)- As consistency 
conditions, we consider Eq. (II3 and [d^l 


V ^PP /pc,pp=0 


1 PcOcp(O), 


(18) 


which involves polymer-colloid correlations. In Figs. IHl and (TUI we also report the bridge functions for 
this case (they are labelled RY-2). For g = 1 we observe a significant improvement with respect to 
the one-parameter RY case, although significant differences with Monte Carlo data are still present for 
fiRg ^ 1- For q — 0.5 instead, the two different RY closures yield equivalent estimates. 

As a final case, we consider the RHNC closure, which relies on the assumption that the bridge functions 
can be accurately parametrized by those of a binary additive hard-sphere mixture. To verify if this is the 
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FIG. 11: Bridge functions for <l?c = 0.3 and q = 0.5 as a fnnction of 6 = rjRg. We report the zero-density function 
obtained by Monte Carlo simulations (MC), and the RY fnnctions for different values of the polymer volume 
fraction <l>p. 


case, we consider q = 0.5 and = 0.3, and compute 

A(i?p) = J \b^p^{r) - b^jf{r,Rp)\ r^dr, (19) 

for different values of the effective polymer radius Rp. The optimal value (minimal A) is obtained for 
Rp = 0.84:2Rg. We can compare this result with that obtained by using the Lado criterion 1 ^. For 
= 0, Eq. ([5]) is satisfied as we use the very accurate hard-sphere correlation function of Ref. [^. To 
determine Rp one needs to consider the linear term in the polymer density, i.e., the equation 

J r^ihgpir) - = 0 . ( 20 ) 

Alternatively, one can determine Rp for several small values of $p, performing at the end an extrapolation 
to ^p —>■ 0. The first method gives Rp = 0.837i?g, while the second one gives Rp = 0.828Rg. Both results 
are very close to the estimate Rp = 0.842i?g obtained by a direct matching of the bridge functions. This 
confirms that the Lado criterion provides the bridge functions that are the best approximations of the 
exact ones. The resulting bridge functions are reported in Fig. (TUI The RHNC estimate of bcp{r) is in 
agreement with the Monte Carlo function for r > 2Rg. As for bpp{r), the RHNC estimate agrees with the 
Monte Carlo one for r > Rg. At smaller distances, instead, the RHNC bridge function underestimates 
the correct one and appears to provide a worse approximation than the RY closure. 

This analysis for <I>p = 0 further confirms the results obtained in Sec. IHI Bl For <I>c = 0.1, the bridge 
functions are quantitatively small, confirming the accuracy of the HNC approximation. On the other 
hand, for <I>c = 0.3, the RY closure is the one that provides the best approximation, while the HNC 
closure is the less accurate one as it cannot reproduce the small-distance behavior of the bridge functions. 
Note that, while bcp{r) is correctly reproduced in the relevant region r ^ Rc, the polymer-polymer bridge 
function is always poorly reproduced for r < Rg. This discrepancy gives rise to similar discrepancies in 
the correlation functions, as discussed in Sec. 11111 


D. Integral equations with Monte Carlo bridge functions 

As a final test we decided to determine the solutions of the integral equations by using the zero-density 
Monte Carlo bridge functions computed in Sec. IHI Cl In other words, we consider the closure relation Q, 
setting for all values of $p, bpp{r; $p) = b^p^(r\ 4>p = 0), bcp{r\ 4>c, 4>j,) = b^'^{r; = 0), and 

brAT\ $c, ‘hp) = <i>c), where the last quantity is the bridge function of a pure hard-sphere system 

|43l |. This approximation is exact for $p = 0 and one may wonder whether it provides a reasonable 
approximation also for $p > 0. We have tested the approach for q = 0.5 and = 0.3. The results 
for the structure factors, reported in Fig. [5] (they are labelled MC-B), show that this approach is only 
marginally better than that based on the HNC closure. Also the termination point, <I>p = 0.11, is only 
slighly above the HNC one, $p = 0.090. 
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To clarify the origin of the discrepancies, we have determined the RY bridge functions for several values 
of <i)p. As the RY estimates reasonably agree with the Monte Carlo data up to the termination line, we 
take them as estimates of the exact density-dependent hap{r\^ci^p)- As one can see from the results 
shown in Fig. [Til the density dependence of the bridge functions is not large (for bcp(r) the relevant region 
is 6 = r/Rg > 2). Yet, this relatively small difference is the cause of the different results obtained. In 
practice, this simple exercise shows that results are extremely sensitive to the specific form of the bridge 
functions in the colloid-liquid phase ^ 0.25. Hence, accurate results can only be obtained by using 
accurate bridge functions, that none of the methods we investigated is able to provide. 


IV. CONCLUSIONS 

In the last years there has been a widespread interest in soft-matter systems characterized by the 
presence of macromolecules of mesoscopic size. In many situations, if one is only interested in the 
thermodynamic behavior or in structural properties on scales much larger than atomic distances, one 
can use coarse-grained (CG) models in which each macromolecule is represented by a single effective 
particle [lainiiii- At variance with simple fluids for which potentials always have a hard core, in 
CG models potentials may be soft, allowing different effective molecules to overlap with a little energy 
penalty. Monocomponent CG models have been extensively studied [E[ni by a variety of techniques. 
Among them, integral-equation methods have been proved to be very accurate. In particular, because of 
the soft nature of the interactions, the HNC and RY closures work quite well [2^, |^. It is then natural 
to investigate whether integral equations can be successfully applied to the study of the phase diagram 
and thermodynamics of more complex systems, for instance to mixtures of macromolecules and colloids, 
characterized by the simultaneous presence of soft and hard-core potentials. 

In this paper, we consider a particular CG model, appropriate to describe long linear polymers inter¬ 
acting with hard-sphere colloids under good-solvent conditionSjU well-studied paradigmatic model whose 
phase behavior has been extensively studied, see, e.g.. Refs, [l^, H^. However, the conclusions should 
have general validity, applying to generic systems with soft and hard-core potentials. The phase diagram 
of the CG model has been discussed recently in Ref. m. The binodal curves and the critical points 
were determined for q = 0.5 and q = 0.8, while, somewhat surprisingly, no sign of phase separation was 
found for g = I up to relatively large polymer densities. Here, we have compared the Monte Carlo results 
with predictions obtained by using integral-equation methods and a variety of different closures: HNC, 
HNC/PY, RY, and RHNC. 

For small values of <i)c we find that HNC is quite succesfull in predicting the correct thermodynamics 
and structure. On the other hand, for = 0.3 (note that the critical point of the fluid-fluid transition 
is located at $c,crit = 0.25 for both q = 0.5 and 0.8) integral equations fail to converge well below the 
binodal line determined by Monte Carlo simulations. Below the termination line the RY closure is the 
one that fares best, reasonably reproducing the zero-momentum structure factors and the pair correlation 
functions. Nonetheless, RY integral equations stop converging at 4)^ = 0.18, 0.39 for <i)c = 0.3 and q = 0.5, 
0.8, respectively, while the binodal is located at significantly larger polymer densities, at 4)^ = 0.38, 0.75 
for the same values of q. 

The failure of integral equations to provide accurate estimates of the phase diagram is probably related 
to the strong nonadditivity of the model. Indeed, similarly large differences are observed in Ref. [l^ for 
systems of nonadditive hard-sphere mixtures. If the system is asymmetric, i.e., for y ;< 0.6 (y is the ratio 
of the diameters of the two spheres, a quantity which is the analog of y), integral equations (and also 
density functional theory) are unable to provide quantitatively reliable results for the phase diagram. 
Moreover, discrepancies increase with the amount of asymmetry considered. 

G.D. acknowledges support from the Italian Ministry of Education Grant PRIN 20I0HXAW77. Gom- 
putations were performed at the Pisa INFN Computer Center and at CINECA (ISCRA PHCOPY 
HPI0CFFG8Q project). 


Appendix A: Technical details 

In the integral-equation approach, pair and direct correlation functions are discretized on N regularly 
spaced points, = nAr. Moreover, all functions are assumed to be zero at a cut-off distance i?niax = 
NAr. Typically, we take Ar = O.OOIRg and N = 32768 or 65536. The grid is extremely hne and 
reasonably large, to guarantee that results are stable with respect to the parameters Ar and N. 
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TABLE II: Estimates of the structure factors Sap{k = 0), of the concentration factor Sc{k), of the virial pressure 
and of the compressibility kt computed using Eq. (1121) for <l>c = 0.3, "Lp = 0.09, q = 0.5, and for the HNC 
closure. We report results for several values of N and Ar. 


Ar = 

N 

= 32768 



N = 

65536 


0.001 

0.002 

0.004 

0.0005 

0.001 

0.002 

0.004 

^p(vir)^3 

0.956 

0.955 

0.953 

0.956 

0.956 

0.955 

0.953 

I3RI/kt 

1.763 

1.762 

1.756 

1.764 

1.763 

1.761 

1.756 

Spp{0) 

3.399 

3.436 

3.560 

3.384 

3.399 

3.436 

3.560 

5cp(0) 

-0.792 - 

-0.801 - 

-0.830 

-0.788 

-0.792 

-0.801 

-0.830 

5cc(0) 

0.263 

0.264 

0.271 

0.261 

0.262 

0.264 

0.271 

5c (0) 

0.396 

0.400 

0.414 

0.393 

0.396 

0.400 

0.414 


In Table mi we report several thermodynamic quantities as a function of Ar and N for the HNC closure 
at = 0.3, <i>p = 0.09, a state point very close to the termination line. 

Estimates do not change as N changes indicating that the cut-off distance is large enough. The step 
size is more crucial, but Ar = 0.001 should be accurate enough. In the paper, most of the analysis use 
Ar = 0.001 and N = 32768. In a few cases, we have checked the results, by changing Ar and/or N by 
a factor of 2. The independence of the results on the chosen parameters allows us to exclude that the 
observed behavior is due either to a too small cut-off distance or to a too coarse discretization of the 
correlation functions. 


Appendix B: Pair potentials 

In this section we report the explicit expressions of the pair potentials. The model consists of coarse¬ 
grained polymers, represented as soft particles, and colloids. Polymers interact via a pair potential Vpp{b) 
given by 

3 

l3Vpp{b) = ^a*exp(-6Vci), (Bl) 

i=l 

where b = rfRg, oi = 0.999225, 02 = 1.1574, 03 = —0.38505, ci = 1.24051, C 2 = 0.85647, and C 3 = 
0.551876. Colloids interact as hard spheres: 

Vcc(.r) = 0 r > 2Rc 

Vcc{r) = +00 r < 2Rc. (B2) 

The polymer-colloid pair potential depends on q. For small values of b = r/Rg, i.e. for b < femin 
(^min ~ Rc/Rg = I/?): the potential pVcp{r; q) is large, hence it is impossible (and practically irrelevant) 
to estimate it accurately. For b > we parametrize it as 

pV,p{r-q) = (B3) 


TABLE III: Coefficients parametrizing fiVcp{r;q) for different values of q. The parametrization is accurate for 
1.91 < b < 5.38, 0.90 < b < 4.54, and 0.47 < b < 4.28 for q = 0.5, 0.8,1, respectively. 


q ai 

ei 

Cl 

02 

62 

C2 

d2 

0.5 0.634486 

0.305183 

2.13936 

15.1368 

0.512611 

1.629090 

1.30679 

0.8 0.411558 

0.318504 

1.40563 

13.5385 

0.728577 0.572266 

1.56655 

1.0 0.982437 

0.496784 

0.98100 

14.1753 

0.84914 

0 

1.6023262 
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where b = r/Rg. Estimates of the coefficients are reported in Table Hill To verify the accuracy of the 
parametrization, we compare the estimate of ^ 2 ,cp = ^ 2 ,cp/^g (^ 2 ,cp is the second polymer-colloid virial 
coefficient) obtained by using the parametrized potential and the estimate of the same quantity in the 
full-monomer model [^. Using the parametrized potentials we obtain A 2 ^cp = 106.79,41.52,27.50 for 
q — 0.5, 0.8, and 1, respectively, to be compared with the full-monomer results ^ 2 ,cp = 107.4(3), 41.7(1), 
27.54(6). Differences are small (they are less than 0.6%), confirming the accuracy of the results. 
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